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O ' 

q | We present a new module of micrOMEGAs devoted to the computation of indirect 

signals from dark matter annihilation in any new model with a stable weakly inter- 
acting particle. The code provides the mass spectrum, cross-sections, relic density 
and exotic fluxes of gamma rays, positrons and antiprotons. The propagation of 
charged particles in the Galactic halo is handled with a new module that allows to 
■ easily modify the propagation parameters. 

1 Introduction 

^ ' Cosmological observations show strong evidence that our Universe contains a large amount 
j> ■ of dark matter (DM). New weakly interacting massive particles (WIMP), such as those 
present in extensions of the Standard Model, have roughly the correct annihilation proper- 
^ ties to fit the high precision cosmological measurements. Several astroparticle experiments 

are actively searching directly or indirectly for this new particle. Indirect detection of dark 
matter particles involves observation of the products of the DM annihilation in the galac- 
tic center, galactic halos or the extra galactic region. The annihilation products include 
positrons, anti-protons, anti-deuterons, gamma-rays and neutrinos. Recently many new 
results from indirect DM searches have been released. Hints of excesses that might be due 
to annihilation of dark matter particles have been reported although an interpretation of 
the measurements in terms of either DM annihilation or some astrophysical source has 
not been confirmed. PAMELA shows an excess in the positron fraction between 10 and 
lOOGeV PP in agreement with earlier indications by HEAT [2] and AMS01 [3J. On the 
other hand PAMELA sees no excess in the antiproton spectrum [I]. Both Fermi [5] 
and ATIC [B] report an excess in the total electron plus positron spectrum but at ener- 
gies of several hundred GeV's, much above those of PAMELA. Furthermore the electron 
spectrum measured by HESS [7] at very high energies is consistent with both Fermi and 
ATIC. The cosmic gamma-rays from the galactic center or from galactic sources have 
been probed in a wide energy range by INTEGRAL [8], Veritas [9], EGRET [TU], as 
well as HESS [ITJ [T2] and Fermi [TBJ. These observations lead to upper bounds on the 
DM annihilation cross section that are however strongly dependent on the halo profile, the 
propagation parameters and the background that is assumed for the standard astrophysics 
processes. Observations in all channels are being pursued actively with in particular Fermi 
and HESS taking data as well as AMS02 to be launched in 2011. 
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The interpretation of the recent and upcoming data requires tools to compute ac- 
curately the signals of DM annihilation in various channels and this in the context of 
different particle physics models. The purpose of the package presented here is to com- 
pute indirect signals in 7, e + and p produced in DM annihilations in the Galaxy. This 
package is presented as a new module of micrOMEGAs [TU [TS], a code that computes the 
dark matter relic density, the elastic scattering cross sections of WIMPs on nuclei relevant 
for direct detection as well as the cross sections and decay properties of new particles rel- 
evant for collider studies]]] micrOMEGAs includes several models of particle physics that 
predict a new stable weakly interacting neutral particle, the minimal supersymmetric 
standard model(MSSM) and several of its extensions, models with extra dimensions or 
the little Higgs model as well as facilities to incorporate new models [IB]. As in earlier 
versions, micrOMEGAs provides the cross sections for dark matter annihilation into SM 
particles and the spectrum for 7, e + and p at the source. The propagation of charged 
particles through the Galaxy which strongly distorts the charged particles spectra is the 
main addition in this version. 

The main features included in the indirect detection module are: 

• Annihilation cross sections for all 2-body tree-level processes for all models. 

• Annihilation cross sections including radiative emission of a photon for all models. 

• Annihilation cross sections into polarised gauge bosons. 

• Annihilation cross sections for the loop induced processes 77 and 7Z in the MSSM. 

• Modelling of the DM halo with a general parameterization and with the possibility 
of taking into account DM clumps. 

• Integrals along lines of sight for 7-ray signals. 

• Computation of the propagation of charged particles through the Galaxy, including 
the possibility to modify the propagation parameters. 

• Effect of solar modulation on the charged particle spectrum. 

• Model independent predictions of the indirect detection signals. 

The neutrino spectrum originating from dark matter annihilation is also computed, 
however the neutrino signal is usually dominated by neutrinos coming from DM capture 
in the Sun or the Earth. The inclusion of this signature is left for a further upgrade. 

In this paper we first review the procedure to obtain the flux of photons or anti- 
particles. This includes the computation of DM annihilation into SM particles and a 
description of the dark matter halo models. The propagation of charged particles includ- 
ing the issue of solar modulation is described in Section 4. The functions available in 
micrOMEGAs are described in section 5 @. Sample results as well as comparisons with 
other codes are presented in Section 6. 

1 DarkSUSY is another public code that computes the indirect signatures of dark matter annihila- 
tion [IT]. It is confined to the minimal supersymmetric model. 

2 micrOMEGAs2.4 contains all the routines available in previous versions although the the format used 
to call some routines has been changed. In particular global variables are used to specify the input 
parameters of various routines. All functions of micrOMEGAs are described in the manual, manual4.tex, 
to be found in the main directory. 
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2 Fluxes from DM annihilation 



Should primordial self-annihilation take place in the early Universe, the same process 
would take place nowadays in the denser regions of the Galactic DM halo. DM annihilation 
in the Galactic halo produces pairs of Standard Model particles that hadronize and decay 
into stable particles These particles then evolve freely in the interstellar medium. The 
final states with 7, e + and p are particularly interesting as they are the subject of indirect 
searches. The production rate of particles from DM annihilation at location x reads 



where av is the annihilation cross-section times the relative velocity of incoming DM 
particles which we evaluate in the limit v = (this is a good aproximation since v = 10~ 3 ). 
Note that (av) includes averaging over incoming particles/antiparticles. For a Dirac 
particel where a xx = 0, (erf) = l/2a xx . m x is the mass of the DM candidate, p(x) is 
the DM mass density at the location x and f a {E) = dN a /dE is the energy distribution 
of the particle a produced in one reaction. The predictions for the energy spectra can 
also depend on some non perturbative effects including QCD and imply the use of Monte 
Carlo simulations. 

2.1 Annihilation cross-sections and energy spectra 

The cross-sections for the different 2-body annihilation channels of WIMPs are calculated 
automatically in any model implemented in micrOMEGAs [16J. This is done through the 
interface with CalcHEP [18J. All cross-sections are computed for a relative velocity v = 0. 
There is an option in the code to define another value for v. The continuum spectrum for 
7, e + , p, v production is calculated as follows. 

The self annihilation of DM particles can occur at tree-level through 16 possible final 
states involving only pairs of SM particles. This includes all flavour diagonal pairs of 
fermions as well as gauge bosons, XX ~~ >* <l<lJ~l + , Wi, W + W~ , ZZ. Other final states 
involving R-parity even particles are also included, in particular the Higgs. In the MSSM 
this corresponds to all the channels with SUSY Higgs particles, namely Zh iy hihj, W ± H Zf 
and H + H~ where hi stands for h, H, A. For the basic channels, qq, W + W~, ZZ and 
gg, we provide tables for 7, e + , p, v production as obtained with PYTHIA version 6.4 [T9] . 
The database has been processed with 2 x 10 6 events at 18 fixed energies corresponding 
to 10 < m x < 5000 GeV. Note that neutrons and antineutrons do not decay in PYTHIA, 
basically a ii is considered to be a p and the small amount of energy lost in the /3 decay 
of ii is neglected. For channels containing two different particles, AB, we obtain the final 
spectrum by taking half the sum of the AA and BB spectra. For channels with Higgses, 
or other particles whose mass are a priori unknown, we recursively calculate all 1 — > 2 
decay channels until we obtain particles in the basic channels. If during these decays we 
get a pair of particles AB where A is one of the basic channel, we suppose that A gives 
half of the spectrum obtained from AA and we continue to decay B. 

The 7 and e + spectra can be substantially modified. First, polarisation of the gauge 
bosons final state can distort the positron and photon spectrum. Second, higher-order 
processes can also significantly modify the particle spectra. For example, photon radiation 
can strongly enhance some channels, this is particularly important for annihilation of a 
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Majorana DM candidate into light fermions which suffers a s-wave suppression. The 
additional photon removes this suppression and the cross section increases by several 
orders of magnitude [20]. The implementation of these two effects in micrOMEGAs are 
discussed in the following subsections. Finally at the one loop level, other final states 
are possible, such as 77, 7Z , gg and 7/1°. Although supressed by a factor a 2 , the 
processes with 7's are nevertheless interesting since they lead to a spectacular signal, a 
mono-energetic 7 ray line. The one-loop processes XX ~~ >" 11 1 1%° have been computed 
automatically with Sloops for the MSSM [2TJ [22], and are incorporated in micrOMEGAs. 
Note that the Majorana nature of the neutralino forbids an annihilation into 7/1° at rest. 

2.1.1 Vector boson polarisation 

The primary particles produced in DM annihilation are by default assumed to be un- 
polarised. In general however these particles and in particular vector particles can be 
polarised. For example the annihilation of neutralinos in the MSSM produces only trans- 
verse W's and Z's while the polarisation of spinor particles can be neglected because of the 
CP invariance of the initial state. The spectra after decay and hadronisation of standard 
particles extracted from PYTHIA also assumes unpolorized particles. 

To take into account the polarisation, we include an option for gauge bosons pair 
production. The first step is to determine the degree of polarisation of the vector bosons 
produced via dark matter anihilation in a given model. We only need to determine the 
polarisation of one of the vector bosons as only the VtVt or VlYl combinations are 
possible. To do this automatically, we compute with micrOMEGAs the three-body process 
XX ~~ >" W~e + v e keeping only the contribution from W pair production followed with the 
W + leptonic decay. We then check numerically the angular distribution of v e in the rest 
frame of the on-shell W + with respect to the direction of flight of the W pair. The angular 
distributions are expected to be 3/8(1 + cos 2 9)dcos9 for Wt and 3/4(1 — cos 2 8)dcos9 
for Wl- With this method, we can reconstruct automatically the W + polarisation in a 
generic model. 

We then need the spectra of the stable particles produced after decay and hadroniza- 
tion of a polarised gauge boson. For this, two methods have been used and compared 
showing perfect agreement. In the first method, we pass to PYTHIA events with four 
outgoing particles representing the decays W + — >■ 2x and W~ — > 2x where the decay 
products are distributed according to the formulas presented above. For this we used 
PYTHIA 6.4 and the Les Houches event interface [23] • Initial events were generated only 
for M x = M\y- Results for heavy CDM where obtained by boosting. The same procedure 
is also used for neutral gauge bosons. In the second method, a reweighting technique is 
applied within PYTHIA 6.4, by measuring event by event the 9 angular distribution of 
the primary W (or Z) decay fermions in the boson rest frame. Namely, for each primary 
fermion a weight is determined depending on the polarisation assumption. For longitudi- 
nal polarisation it is equal to | x (1 — cos 2 9) while it becomes | x (1 + cos 2 9) for transverse 
polarisation as explained previously. Then for each stable particle (7, e + p, u), the weight 
of the W ( or Z) decay fermion they originate from, is used to build the energy spectrum 
for each polarisation scenario. This is done for 18 different CDM masses. Extensive com- 
parisons of spectrum distributions for direct or reweighted longitudinally polarized bosons 
have been made. They are in perfect agreement. We also checked that the average of the 
polarised spectra obtained with each method agreed with the unpolarised one. 



4 



Taking into account the polarisation of the W's (or Z's) leads to a harder spectrum 
for e + originating from transverse W's than if the W's were assumed unpolarised. The 
polarisation effect corresponds to at most a factor of 3/2 increase for the most energetic 
charged particles. This is because both transversely polarised W give a harder spectrum 
while the spectrum for longitudinally polarised W vanishes at high energy, see fig. [T] for 
the positron spectrum. Both the polarised and unpolarised spectra are available in the 
basicSpectra routine. To calculate the spectra of DM annilation taking into account 
the polarisation of vector bosons one has to set a switch in the calcSpectrum routine, 
see the routines description in Section HI 
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Figure 1: dN/dE for positrons from DM annihilating into W + W~ with transversely 
polarised (full/black), longitudinally polarised (dash/pink) or unpolarised W's (full/blue), 
here m x = 1 TeV. 

2.1.2 Photon radiation 

In a DM annihilation process, photons can be emitted from the external legs, final state 
radiation (FSR), from an internal leg, from a quartic vertex or in the subsequent decay 
of the pair of particles directly produced. The first two processes cannot be treated 
separately in a gauge invariant manner and are associated with a three-body annihilation 
process. They are a priori suppressed by a factor a, however large enhancements can 
occur [21]. In particular final state radiation features a collinear log corresponding to 
a photon emission from an external leg. This virtuality gives a factor log(4M^/mj) 
where mj is the mass of the outgoing particle. For light particles final state radiation 
can therefore be approximated by the production cross section for two particles times 
a radiation factor that only depends on the spin of the particle. This is the approach 
followed in PYTHIA. Furthermore in the case of light masses rrif M x one expects 
multiple photon radiation for the soft spectrum, this is also used in PYTHIA so that the 
FSR spectrum differs from the one corresponding to a single photon radiated from an 
external leg by up to 10% - 20% in the high energy part. For W bosons the virtuality 
is never large so photon radiation is ignored in PYTHIA, the only photons included in 
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PYTHIA are those coming from the decays of the gauge bosons. In this case a full 
2—7-3 calculation is in order, however barring exceptional cases the yield is suppressed 
by a factor a (with no large logarithm enhancement) and therefore small. An important 
example where the approximation, and therefore applying PYTHIA, fails is if the cross 
section without radiation of a photon is for symmetry reasons very small or vanishing so 
that factorisation does not hold. A notable case, mentionned above, is the annihilation 
of the Majorana neutralino, at v = 0, to a pair of almost massless fermions. The s-wave 
cross section of order the chirality factor m 2 /M 2 can be much smaller than the radiation 
cross section of order a since the chirality argument no longer applies once a photon is 
radiated. The cross section including a photon emission can be enhanced in some specific 
cases. This can occur in situations where the t-channel particle exchanged between the 
DM particle and the external charged particle is not far from the DM mass, leading to an 
enhancement factor M 2 /(Mf — M 2 ). Mi is the mass of the internal particle. 

In our code we compute the direct photon production through the full 2 — > 3 cal- 
culation, those are generated at run-time but only for the situation of interest which 
we have defined as M| < 1.5M 2 . We also compute the full WW'j in situations where 
M x > 500 GeV because of the potential large log enhancement due to photon emission 
since in this situation the W can be relativistic. All 3-body final states are included in 
the computation of the photon spectrum when this option is chosen. On the other hand 
only the e + e~7 process is included in the positron spectrum. Indeed these are the only 
processes that affect significantly the hard part of the positron spectrum. 

When the full XX ~>* ffl process is generated by the code, one should avoid double 
counting due to the inclusion of the photons induced by PYTHIA through the direct 
ffj. By direct we mean the photons generated prior to the possible decays of / and 
hadronisation. To subtract the FSR contribution already taken into account in PYTHIA 
we remark that our generated XX ~* ffl photon must exhibit the infrared divergent 
behaviour 1/x, if the XX ~~ >* I f cross section is not vanishing. This infrared behaviour at 
x = can only originate from radiation from external charged legs. We therefore expand 
the generated 2 — ^3 cross section of our code dav/dx around x = and write 

dav A „ _ 

— = - + B + Cx 2 

dx x 

The idea is to subtract from our generated cross section the A/x term obtained for small 
enough x. For this we use x = 0.01. We have checked that 0.001 < x < 0.03 give 
similar results. There might still be a mismatch between the coefficient A extracted from 
the full calculation and the one contained in PYTHIA. This difference is due to QCD 
hadronisation, higher order terms, additional photons and choice of the scale Q 2 for the 
splitting, contained in PYTHIA. To conform with the splitting function, for light fermions 
for example, we coud then subtract A(l — x + x 2 /2)/x, however in our code this hardly 
makes a difference in situations where the full calculation is important. At the level 
of implementation let us mention that there might be a problem caused by the finite 
precision of the phase space integration of the 2 — > 3 matrix elements. For a fixed photon 
energy, there is an infrared pole when the outgoing fermion and the radiated photon 
become collinear. Furthermore strong numerical cancellations occur when summing over 
all diagrams. In the case of small fermion masses this can lead to numerical instability. 
To solve this we replace the integration variable d cos (ft to deft ((ft is the angle between the 
outgoing fermion and the radiated photon in the rest frame of the fermion pair). This 
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trick works well for nif < 10~ 4 m x . When irif < 10~ 3 m x we keep g?cos0 as the integration 
variable and impose a cut | cos<fr\ < 0.99, in this region the numerical calculation is robust. 
The poles at <fi = rrif/M x due to radiation from the final charged particles are safely out 
of the integration region. We have tested that both methods are in good agreement when 
rrif = 10 _3 m x . 

The importance of the extra contributions not contained in the factorised FSR to 
the photon spectra is illustrated in Fig. [2^ for the CMSSM point mo = 70 GeV, Mi = 
250 GeV, A = -300 GeV, tan/3 = 10. Here M x = 97.9 GeV, M f = 107.6 GeV, 
Me = Mfi = 123.9 GeV, neutralinos annihilate dominantly into tau pairs, this is an ex- 
ample of a case where the t-channel particle exchanged, the f, is not far from resonance 
so that photon radiation from internal lines is enhanced. Note that the curve denoted 
TT7 in Fig. [2k, represents only the contribution from the 3-body process, the additional 
photons that originate from r decays are included when computing the full process. The 
enhancement in the high-energy part of the photon spectrum in a case where DM annihi- 
lation into gauge bosons pairs is dominant is illustrated in Fig. [2b for the MSSM. Here the 
DM is a mixed bino/higgsino LSP with fi = 545GeV, M 1 = 500GeV, M 3 = 6M 1 = 3M 2 , 
tan/3 = 20, Ma = 2TeV and all soft sfermion masses heavy, m? = 2.5TeV. 




Figure 2: a) Photon spectrum for the CMSSM point m = 70 GeV, Mi = 250 GeV, 
Aq = —300 GeV, tan /3 = 10 for XX ~^ t + t~7 (black), e + e~7 (black-dash), with photons 
from external legs from PYTHIA (FSR) (green-dash) and all contributions from 3-body 
final states (green-full), b) Photon spectrum for a MSSM point (// = 545 GeV, Mi = 
500GeV,tan/3 = 20, Ma = 2TeV and m~ f = 2.5TeV) with (full) and without (dash) the 
contributions from the three-body process XX ~^ W + W~ , j 



2.2 Dark matter halo models 

In micrOMEGAs routines which calculate the propagation of particles in the Galaxy the 
DM halo distribution is an input parameter. Thus micrOMEGAs can work with any 
sphericaly symmetric DM halo profile. As an example of DM halo distribution we include 
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a widely used spherically symmetric parametrization of the dark matter halo 



p(r) 

Fhaloir) 



P&F halo (r) 

r ]7 [1 + (r Q /a) a 
1 + (r/a) a 



(3) 



The values for the ft, 0, 7 and a parameters for the most common halo models are listed 
in Tab. [TJ the default values are those of the NFW profile. The value of the DM density 
at the solar location p Q and r Q the distance of the Sun to the Galactic center are global 
parameters of micrOMEGAs. Their value can be redefined simply and the default values 
are listed in Table 3. An important remark concerns the central divergence. For 7 > 1.5 
in eq. [3] there is an non integrable squared density in the center of the galaxy. To avoid 
this singularity we set a limit for the distance from the Galactic center r > r min = O.OOlpc 
in our integration routines. This is done for any halo profile. 



Halo model 


a 





7 


a (kpc) 


Isothermal with core 


2 


2 





4 


NFW 


1 


3 


1 


20 


Moore 


1.5 


3 


1.5 


28 



Table 1: Halo parameters for three common profiles 

The user also has the possibility to use a totally different parameterization for the halo, 
provided it is spherically symmetric. For example the Einasto profile has been advocated 
recently [25] 

-2 ( ( r 



Fhaloir) = exp 



— 11 — 1 -1 



(4) 



where the default value for a is 0.17. 

DM annihilation depends on the squared density p 2 {r). In general a clumpy structure 
of DM will lead to p 2 (r) > p(r) 2 . The effect of clumping can be described by another 
profile 
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The presence of clumps will lead to an enhancement factor or boost factor. In general the 
clump number density in a volume bounded by the characteristic diffusion length of the 
involved species will determine the size of the enhancement factor. The resulting boost 
factor can therefore differ for photons, positrons or antiprotons. In the last two cases it 
does not typically exceed 20 [26] . Larger boost factors can be found for some extreme 
clump configurations, for example a big sub-halo close to the Earth [26], although this 
situation is very unlikely 



2.3 Photons 

The gamma ray flux can be evaluated as 
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m 
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and is expressed in number of photons per cm 2 s sr. The factor H includes the integral 
of the squared of the dark matter density over the line of sight, 

where r' = ^/r 2 + r| — 2rr Q cos and is the angle in the direction of observation. In 
micrOMEGAs one can compute the photon flux performing the integral over the line of 
sight and over the opening angle which characterizes the detector resolution, see section 21 

3 Galactic propagation of charged particles 
3.1 General framework 

The charged particles generated from DM annihilation propagate through the Galactic 
halo and their energy spectrum at the Earth differs from the one produced at the source. 
Charged particles are deflected by the irregularities of the galactic magnetic field. In 
the Milky Way which has strong magnetic turbulence Monte Carlo simulations indicate 
that this is described by space diffusion. Charged particles also suffer energy losses from 
synchroton radiation and inverse Compton scattering as well as diffusive reacceleration in 
the disk. Finally galactic convection wipes away charged particles from the disk. Solar 
modulation can also affect the low energy part of the spectrum. The equation that 
describes the evolution of the energy distribution for all particles (protons, anti-protons, 
positrons) reads 

g- z {VciJa) - V • (K(E)V1> a ) - (KE)^a) = Qa(x, E) (8) 

where ip a = dn/dE is the number density of particles per unit volume and energy, a 
denotes the particle specie and Q a is the source term. Here we do not consider back- 
ground and include only particles produced from dark matter annihilation, Eq. [TJ For 
antiprotons a negative contribution to the source term will also be considered to account 
for antiproton annihilations in the interstellar medium (see section 1X5]) . K is the space 
diffusion coefficient, assumed homogeneous. 

K(E) = K f3(E) (K/l GV) S (9) 

where /3 is the particle velocity and 1Z — p/q its rigidity. b(E) is the energy loss rate. 
Another term that describes the diffusive reacceleration has been neglected. The simple 
power law for K(E) is inferred from magnetohydrodynamics considerations [29], once 
the convective velocity is taken into account it has been shown that this form for K(E) 
was adequate to fit the B/C data [30J. 

The propagation equation, Eq. [HI is solved within a semi-analytical two-zone model 
which was discussed in [28[ [30], [32]. Within this approach the region of diffusion of 
cosmic rays is represented by a thick disk of thickness 2L and radius R ~ 20 kpc. The 
thin galactic disk lies in the middle and has a thickness 2h ~ 200 pc and radius R. 
The boundary conditions are such that the number density vanishes at z = ±L and at 
r = R. The galactic wind is directed outward along the z direction so the convective 



drp 2 (r') 



(7) 
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velocity is also vertical and of constant magnitude Vc(z) = Vcsign(z). The propagation 
parameters 5, Kq, L, Vc are constrained by the analysis of the boron to carbon ratio, a 
quantity sensitive to cosmic ray transport [30]. Typical values for these coefficients are 
listed in Table El default values are those of the MED model. 



Model 


8 


K (kpc 2 /Myr) 


L (kpc) 


V c {km/s) 


MIN 


0.85 


0.0016 


1 


13.5 


MED 


0.7 


0.0112 


4 


12 


MAX 


0.46 


0.0765 


15 


5 



Table 2: Typical diffusion parameters that are compatible with the B/C analysis [30| I3T] 

The solution for the energy distribution, eq. [HI will generally be expressed as an integral 
equation 



dE s / ^x s G(x ,£;x 5 ,£ 5 )Q(x 5 ,£c 



(10) 



where the last integral is performed over the diffusive halo. Eg is the energy at the source 
and G(x , E; xg, Eg) is the Green function which determines the probability for a cosmic 
ray produced at X5 with energy Eg to reach a detector at the Earth with an energy E. 
The differential flux is related to the number density 



^ = IT*- 



'in 



where v is the particle velocity. The method of solution of the general equation adapted 
to each type of charged cosmic rays will be described next. 



3.2 Positrons 

The energy spectrum of positrons is obtained by solving the diffusion-loss equation keeping 
only the two dominant contributions: space diffusion and energy losses. 

- V • (K(E)VM - A (b{E)i, e+ ) = Q e +(x, E) (12) 

Here K = K (E / E ) 5 since for energies above 0.1 GeV the positrons are ultra relativistic 
and the rigidity 1Z is proportionnal to E, E Q = lGeV. The positron loss rate is dominated 
by synchrothon radiation in the galactic magnetic field and inverse Compton scattering 
on stellar light and CMB photons, 

KE) = (13) 

where te = 10 16 s is the typical energy loss time. Subdominant contributions from 
diffusive reacceleration and galactic convection were shown to have an impact only below 
a few GeV's [38], these effects are not included in our treatment. 

The propagation equation can be transformed into the heat conductivity equation 
after the substitutions 

tl>{E,r,z) = N(t,r,z)/b(E) (14) 
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The propagation equation now reads 
(|-V W ,r,,) = 



(16) 



E=E(t) 



Assuming p 2 (r, z) = p 2 (r,—z), we can solve Eq. [TH] in the region < z < L with the 
boundary conditions 

dN(t,r,z) 



N(t,r,L) = 0; 



dz 







(17) 



z=0 



The Green function for Eq. [THlcan be factored into a r-dependent part, a Gauss function, 
and a ^-dependent part. The latter is more complicated because of boundary conditions, 



J^-V 2 )G(r,r,z,z') = 5(r)5 2 (r}5(z - z') 



G(t, r,z,rf) 



©(r) ( -r\ 



d d 2 



(18) 
(19) 
(20) 



To calculate g v , the vertical Green function, we construct a basis of eigenstates e n (z 

dz 2 



of the operator — defined on an interval < z < L with the boundary conditions 



e n (L) = and e' n (0) = 0. This basis is 

e n (z) = sin(7r(n + 0.5)/L) 
The vertical Green function expressed in terms of e n (z) reads 

oo 

g v (t,z,z') = Q(t) ^2e- k ' t e n (z)e n (z')/c n 



(21) 



(22) 



n=0 



where c n are the normalization constants 

-L 



(-7i\Z) (-m. \Z) dz 5 nm c n 



(23) 







When t is small, eq. [22] does not converge well. In this case g v is rather obtained using 
the method of electrical images, 



II — — r*n \ 



1)™ / (z-z'+2nL) 2 {z + z' + 2nL) 2 

e 4t + e 4t 



(24) 



Using the Green functions, Eq. [T8]-[20], the positron density near the Sun is obtained after 
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a 3D integration 



MEo,r Q ,0) = J dEf(E)D(t(E )-t(E),r Q ) (25) 

Eo 



1 



D(0,r Q ) = — p2(r ,O) (26) 

L oo 

D(r,r e ) = l/^(x, z ,0)/rdr^exp(-t^)!)/,(M) (27) 







Ji(x) = - / d(f)e- 2xsin2 2 =I (x)e- x (21 

where Jo is the modified Bessel function of the first kind. In practice we use R as the upper 
limit of the integral over r in Eq. [271 A more precise treatment of the boundary condition 
is not required as the positrons originating from far away sources suffer significant energy 
losses. 

Note that D(r, r ) is an universal function for all energies. To calculate il>g(E) for all 
positron energies, it is more efficient to first tabulate D as a function of r in the region 
< r < t(E m i n ) — t(M x ) and then perform a fast integration for all energies. This is the 
method implemented in the micrOMEGAs routine posiFluxTab. 



3.3 Antiprotons propagation 



The propagation of antiprotons is dominated by diffusion and the effect of the galactic 
wind. The source term includes the annihilation of DM, Eq. [1], as well as a negative source 
term corresponding to the annihilation of antiprotons in the interstellar medium (H, He). 
The annihilation rate 

T tot = a; n H n v p n H + a^eVpnne (29) 

where Vp is the velocity of the p. The annihilation cross-sections cr^ n are found in [Ml 135] 
and rescaled by a factor 4 2 / 3 for <r|^. The average densities in the galactic disc are set 
to fin = 0.9cm -3 and = 0.1cm -3 . The production of secondary antiprotons is not 
included in the code. 

The energy spectrum of antiprotons is obtained by solving the diffusion equation 



-K(E)V 2 + V c ^ + 2(V C + hT tot (E))5(z) 



ipp(E,r,z) 



<tv p 2 (r, z) 



Ml 



fp{E) (30) 



An important difference with the positron case is that energy loss of antiprotons is negli- 
gible, we will therefore omit the E dependence until the final formula. dip(r, z)/dz should 
have a discontinuity at z = 0. Assuming that p 2 (r, z) = p 2 (r, —z) means that ipp(r, z) has 
the same symmetry and the discontinuity can be presented as a boundary condition 



djj p (r, z) 
dz 



(31) 



2 = 
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After substituting 

ipp(r, z) = exp(k c z)N(r, z); (32) 
where k c = V C /(2K), Eq. [301 simplifies to 



[-V 2 + k 2 ] N(r, z) = e-^^ P ^f P (33) 



for < z < L with the boundary conditions 

dN(r, z) 



N(r,L) = 0; 



dz 



= N(r,0)v d - (34) 

z=0 



where v d = (V c /2 + hT tot ) / K. 

To compute the Green function for Eq. f)33p we proceed as for the positron case and 
construct a basis of eigenstates of the — operator defined on the interval < z < L 
with the boundary conditions 

e n (z) = sm(k n (L - z)) (35) 
where the set of k n is defined by the condition 

<(0) = e n (0)v d . (36) 

The Green function then reads 

1 oo 

G(r,z,z') = Y, K ^^W+K)zn{z)e n {z')/c n (37) 



n=0 



where c n are the normalization constants (see Eq. 123]) and Kq is the MacDonald function 
defined by the equation 

A r K {r) - K (r) = -2n5 2 (r) (38) 
The antiproton energy spectrum at the Earth is obtained after a 3D integration 



ME, r Q , 0) = rdr jf ^G(r, z, 0)e" fc - £ (39) 



with r' = + r 2 + 2r r cos 0. 

To avoid numerical problems due to a possible singularity in the dark matter density 
near the center of the Galaxy, we integrate Eq. [39] in the central region |x| < tq = O.Olkpc 
fixing p 2 (|x|) = p 2 (r ), we then treat the \x\ < r region as a point-like source with 
p 2 (|x|) — p 2 (r ). The radial boundary conditions can be simulated by modifying the 
DM density. First we note that sources located at a distance r > R have a negligible 
effect, we then take p = when r > R. Second, radiation from a source located near 
the boundary r < R will be suppressed. To estimate this suppression we add mirror 
sources with negative charges on each side of the boundary r = R. The long distance 
contribution from a point-like source is proportionnal to K (r/a) ~ exp(— r/a) where 
a = 1/y/kl + k 2 c w 2L/tt. We then modify the DM density 

p 2 (r,z) e(R-r)(l-exp[-2(R-r)/a]) p 2 (r,z) (40) 
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The suppression factor is significant only in the region R — L < r < R. 

The integrand in Eq. ( 139]) (1(E)) is a smooth function that features only a slight 
energy dependence through the terms K(E) and T to t(E). To provide a fast calculation of 
the antiproton spectrum we interpolate the function 1(E) in the range M x — -Emm, and 
multiply this interpolated function by fp{E) to obtain the final result for all E. E min is 
defined by the user. Note that points for interpolation are added automatically until the 
required precision is reached. 

3.3.1 Solar modulation 

When charged particles propagate through the solar system, they are affected by the solar 
wind and lose energy. This effect leads to a shift in the energy distribution between the 
interstellar spectrum and the spectrum at the Earth, this shift affects the low energy part 
of the spectrum. We implement solar modulation using the force field approximation [55] . 
In this approximation the flux at the Earth is simply related to the flux at the heliospheric 
boundary (<f> h ), 

dE Q V \dE h 1 ' 

where p & and pt are the momenta at the Earth and the heliospheric boundary. The 
energies at the two locations are simply related by 

E Q = E IS - \Z\<f> F , (42) 

where the Fisk potential 4>f varies between 300MV to 1000MV from the minimum to the 
maximum of solar activity. The default value is set to 500MV. 



4 Functions of micromegas 

A complete description of all micrOMEGAs functions is available in the online manual, 
http://lapth.in2p3.fr/micromegas/. This updated manual is also provided with the 
code, the file manual4.tex can be found in the main micrOMEGAs 2.4 directory. Here we 
describe the functions that are specific to the indirect detection module. A new feature of 
version 2.4 is the use of global parameters. A list of the global parameters relevant for the 
indirect detection module and their default values are given in Table 1. The numerical 
value for any of these parameters can be simply reset anywhere in the code. 

4.1 Spectra interpolation and display 

Various spectra of SM particles produced in DM annihilation processes are stored in arrays 
containing NZ=250 elements. The i th element of an array corresponds to dN/dzi where 
Zi = \og(Ei/M x ). Here E_i is the kinetic energy of theparticle and N stands for either 
the number of particles or a particle flux in (cm 2 s sr) _1 o 

The following functions can be used for interpolation and visualization 



3 Although not needed if one uses the interpolation functions, the value of Zi can be obtained by the 
function Zi(i). In the current version Zi(i) = log(10~ r )(i/NZ) 15 . 
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Name 


Default value 


Units 


Comments 


Mcdm 




GeV 


Mass of Dark Matter particle, M x 


Rsun 


8. 


kpc 


Distance from the Sun to the Galactic center, r 


rhoDM 


0.3 


GeV / cm 3 


Dark Matter density at Rsun, p Q 


Rdisk 


20 


kpc 


Radius of the galactic diffusion disk, R 


K_dif 


0.0112 


kpc 2 /Myr 


Diffusion coefficient Kq 


L_dif 


4 


kpc 


Half height of the galactic diffusion zone L 


Delta_dif 


0.7 




Slope of diffusion coefficient, 5 


Tau_dif 


10 16 


s 


Positron energy loss time scale, te 


Vc_dif 


12 


km/s 


Convective velocity of Galactic vind , Vc 



Table 3: Global parameters of the indirect detection module 



• SpectdNdE(E, spectTab) 

interpolates the tabulated spectra and returns the dN/ dE distribution where E is the energy 
in GeV. 

• zlnterp(z, SpectTab) 

interpolates the tabulated spectra and returns the dN/dz distribution where z = log(E/M x ), 
here z = corresponds to E — M x . 

• displaySpectrum (Spectrum, message ,Emin,Emax, Units) 

displays the resulting spectrum, message is a text string which gives a title to the graphic 
plot. Emin and Emax define energy cuts. If Units=0 the spectrum is written as a function 
of z otherwise the spectrum is a function of the energy in GeV. 

4.2 Annihilation spectra 

• calcSpectrum (key , Sg , Se , Sp , Sne , Snm , Snl , feerr) 

calculates the spectra of DM annihilation at rest and returns av in cm 3 / s . The calculated 
spectra for 7, e + , p, v e , z/ M , v T are stored in arrays of dimension NZ as described above: Sg, 
Se, Sp, Sg, Sne, Snm, Snl. To remove the calculation of a given spectra, substitute NULL 
for the corresponding argument, key is a switch to include polarisation of W,Z bosons 
(key=l) or photon radiation (key=2). Photon radiation is added to all subprocesses when 
computing the photon spectrum while only the 3-body process XX ~^ e + e~7 is included 
for the positron spectrum. When key=4 the cross sections for each annihilation channel 
are written on the screen. More than one option can be switched on simultaneously 
by adding the corresponding values for key. For example both the W polarisation and 
photon radiation effects are included if key=3. A problem in the spectrum calculation 
will produce a non zero error code, err 7^ 0. 

• spectrlnf o (Xmin, spectrTab,&Ntot ,&Xtot) 

provides information on the spectra generated. Here Xmin defines the minimum cut for 
the energy fraction x=E/Mcdm, Ntot and Xtot are calculated parameters which give 
on average the total number and the energy fraction of the final particles produced per 
collision. Note that the upper limit is Xtot=2. 

• basicSpectra(pdgN,outN,Spectr) 

is a routine for model independent studies that computes the spectra of outgoing particles 
as obtained by PYTHIA and writes the result in an array of dimension NZ, Spectr. pdgN 
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is the PDG code of the particles produced in the annihilation of a pair of dark matter 
particles. The spectrum for polarised W's or Z's is obtained by substituting pdgN+ , T ) 
(transverse) or pdgN+'L' (longitudinal) for the PDG code of the W(Z), pdgN=24(23). 
outN specifies the outgoing particle, 

outN = {0, 1,2,3,4,5} for {7, e + ,p~, u e , u T } 

Note that the propagation routines for e + ,p, 7 can be used after this routine as usual. 

• loopGamma(&vcs_gz ,&vcs_gg) 

calculates av in cm 3 / s for the loop induced neutralino annihilation into 7Z and 77. In 
case of problem the function returns a non-zero value. This function is available only for 
the MSSM. 

4.3 Distribution of Dark Matter in the Galaxy. 

To compute the signal from an indirect detection experiment one has to take into account 
the dark matter distribution in the Galaxy. Both the DM density profile as well as the 
clump profile, Eq. [31 have to be defined. The DM density at the Sun, p as well as r , 
the distance from the center of the Galaxy to the Sun are defined by the global variables 
rhoDM and Rsun. 

• setHaloProf iles(F halo , F dump ) 

allows to change both the halo and the clump profile. Any sphericaly symmetric DM halo 
profile can be defined. 

• hProf ileABG(r) 

is the default halo density profile, Eq. [3] 

• setProf ileABG ( alpha , beta , gamma , a) 

resets the parameters of the DM distribution profile. The default parameters correspond 
to the NFW profile, a = 1, p = 3, 7 = 1, a = 20[kpc]. 

• hProf ileEinasto(r) 

is the Einasto halo density profile, eq. HJ 
•setProf ileEinasto (alpha) 

sets the parameter a for the Einasto profile, the default value is a = 0.17. 

• noClumps(r) 

is a non clumpy profile which is used by default. It returns 1 for any argument. 

4.4 Particle propagation. 

The spectrum of charged particles observed strongly depends on their propagation in the 
Galactic Halo. The propagation depends on the global parameters 

K_dif , Delta_dif, L_dif, Rsun, Rdisk 

as well as 

Tau_dif (positrons) , Vc_dif (antiprotons) 

• posiFluxTab(Emin, sigmav, Se, Sobs) 

computes the positron flux at the Earth. Here sigmav and Se are values obtained by 
calcSpectrum. Sobs is the positron spectrum after propagation. Emin is the energy cut 
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to be defined by the user. Note that a low value for Emin increases the computation time. 
The format is the same as for the initial spectrum. The function SpectrdNdE(E,Sobs) 
described above can also be used for interpolation, in this case the flux returned in (GeV 
s cm 2 sr) _1 ). 

• pbarFlux(E,dSigmavdE) 

computes the antiproton flux for a given energy E and a differential cross section for 
antiproton production, dSigmavdE. For example, one can substitute 
dSigmavdE=cxt;SpectdNdE (E , SpP) 

where av and SpP are obtained by calcSpectrum. This function does not depend on 
the details of the particle physics model and allows to analyse the dependence on the 
parameters of the propagation model. 

• pbarFluxTab(Emin, sigmav, Sp, Sobs) 

computes the antiproton flux, this function works like posiFluxTab, 

• solarModulation(Phi , mass, stellarTab, earthTab) 

takes into account modification of the interstellar positron/antiproton flux caused by the 
electro-magnetic fields in the solar system. Here Phi is the effective Fisk potential in 
MeV, mass is the particle mass, stellarTab describes the interstellar flux, earthTab is 
the calculated particle flux in the Earth orbit. 

The photon flux does not depend on the diffusion model parameters but on the angle 
between the line of sight and the center of the galaxy as well as on the annihilation 
spectrum into photons 

• gammaFluxTab(f i , df i , sigmav, Sg, Sobs) 

multiplies the annihilation photon spectrum with the integral over the line of sight and 
over the opening angle to give the photon flux, f i is the angle between the line of sight 
and the center of the galaxy, df i is half the cone angle which characterizes the detector 
resolution (the solid angle is 2n{l — cos(dfi)) , sigmav is the annihilation cross section, 
Sg is the DM annihilation spectra. Sobs is the spectra observed. Note that Emin is not 
specified, since this function is not time consuming the integration is done for all energies. 
The function gammaFluxTab can be used for the neutrino spectra as well. 

• gammaFlux(f i ,df i ,vcs) 

is the same function as gammaFluxTab above but corresponds to a discrete photon spec- 
trum, vcs is the annihilation cross section, for instance in the MSSM it is calculated with 
the loopGamma function. The function returns the number of photons per cm 2 of detector 
surface per second. Note that for XX ~~ 77 the result should be multiplied by a factor 2 
as each annihilation leads to the production of two photons. 

Note that for solarModulation and for all *FluxTab routines one can use the same 
array for the spectrum before and after propagation. 

5 Examples and results 
5.1 Sample output 

Running the main.c sample file in the micrOMEGAs/MSSM directory choosing the options 
MASSES. INFO, CONSTRAINTS, OMEGA, INDIRECT.DETECTION, CDM.NUCLEON will lead to the 
following output. 
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========= SLHA file input ========= 

Initial file "model2 . slha" 
Warnings from spectrum calculator: 
Model: model2 

Dark matter candidate is '~ol' with spin=l/2 mass=l . 48E+02 

"ol = 0.833*bino -0.114*wino -0.448*higgsinol -0 . 303*higgsino2 

=== MASSES OF HIGGS AND SUSY PARTICLES: === 
Higgs masses and widths 
h 123.08 2.59E-03 
H 1000.20 1.19E+01 
H3 1000.00 1.20E+01 
H+ 998.63 1.18E+01 



Masses of odd sector Particles: 



~ol 


MNE1 




147 


7 I 


I "2+ 


MC2 




189 


9 1 


I ~o2 


MNE2 




198 


2 


~o3 


MNE3 




211 


1 | 


I ~o4 


MNE4 




345 


1 | 


I "1+ 


MCI 




345 


3 


~g 


MSG 




1108 


9 I 


I ~tl 


MStl 




2418 


8 I 


I ~bl 


MSbl 




2496 


4 


~ne 


MSne 




2499 


2 I 


I ~nm 


MSnm 




2499 


2 I 


1 ~nl 


MSnl 




2499 


2 


~uL 


MSuL 




2499 


4 I 


I ~cL 


MScL 




2499 


4 I 


I ~uR 


MSuR 




2499 


8 


~cR 


MScR 




2499 


8 I 


I ~11 


MS11 




2499 


8 I 


I ~sL 


MSsL 




2500 


1 


~dL 


MSdL 




2500 


1 | 


I ~mL 


MSmL 




2500 


3 I 


I ~eL 


MSeL 




2500 


4 


~eR 


MSeR 




2500 


4 I 


I ~mR 


MSmR 




2500 


5 I 


I ~dR 


MSdR 




2500 


7 


~sR 


MSsR 




2500 


7 I 


I ~12 


MS12 




2501 


I 


I ~b2 


MSb2 




2504 


4 


~t2 


MSt2 




2589 


4 I 























==== Physical Constraints: ===== 

deltartho=5 . 70E-06 

gmuon=-8.31E-ll 

bsgnlo=3.94E-04 

bsmumu=3 . 02E-09 

btaunu=9.96E-01 

MassLimits OK 

==== Calculation of relic density ===== 
Xf=2.45e+01 0mega=1.13e-01 

Channels which contribute to 1/ (omega) more than 1%. 
Relative contrubutions in 7, are displyed 
60% "ol ~ol -> W+ W- 
26% "ol "ol -> Z Z 
97, "ol "ol -> Z h 
17, "ol "ol -> h h 
The rest 3.28 7, 

==== Indirect detection ======= 

Channel vcs[cs"3/s] 
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~ol,~ol -> h h 1.30E-43 
~ol,~ol -> Z h 1.31E-27 
~ol,~ol -> Z Z 5.03E-27 
~ol,~ol -> W+ W- 1.20E-26 
sigmav=l . 83E-26 [cm~3/s] 

Photon flux for angle of sight f=0.00[rad] 
and spherical region described by cone with angle 0.04[rad] 
Photon flux = 2.0 IE- 14 [cm" 2 s GeV] "{-1} for E=73.8[GeV] 
Gamma ray lines : 

E=1.34E+02[GeV] vcs(Z,A)= 1 . 37E-29 [cm~3/s] , f lux=7 . 19E-14 [cnT2 s] "{-1} 
E=1.48E+02[GeV] vcs(A,A)= 2 . 47E-30 [cm~3/s] , f lux=2 . 58E-14 [cm~2 s] "{-1} 
N_=19 

Positron flux = 8 . 09E-12 [cm~2 sr s GeV] "{-1} for E=73.8[GeV] 
N_grid=9 

Antiproton flux = 1 . 77E-11 [cm"2 sr s GeV] "{-1} for E=73.8[GeV] 

==== Calculation of CDM-nucleons amplitudes ===== 
CDM-nucleon micrOMEGAs amplitudes: 
proton: SI -4.570E-09 SD -6.116E-07 
neutron: SI -4.582E-09 SD 5.355E-07 
CDM-nucleon cross sections [pb] : 

proton SI 9.015E-09 SD 4.844E-04 

neutron SI 9.061E-09 SD 3.714E-04 

5.2 Comparison with other packages 

We have performed many tests to check the consistency of the results obtained with 
micrOMEGAs. For the treatment of the hadronization process, we have computed the 
fragmentation functions with PYTHIA(those are implemented in micrOMEGAs) and with 
HERWIG. Despite the fact that the two codes are based on different modelling of the 
hadronization process, the functions dN/dE obtained with PYTHIA and HERWIG are very 
similar. A comparison between the two codes in the case of a 500 GeV annihilating into 
bb pairs has been performed and gave similar results. 





Model 1 


Model 2 




H = -440, M A = 1000 


H = -200, M A = 1000 




M 2 = 800, M = 2500 


M 2 = 320, M = 2500 




A t = A b = 1000 


A t = -At = -2500 




tan/3 = 10 


tan /3 — 8. 




micrOMEGAs 


DarkSusy 


micrOMEGAs 


DarkSusy 


M x (GeV) 


386.7 


386.7 


147.7 


147.7 


crv(cm 3 s _1 ) 


1.97 x 10" 26 


2.32 x 10" 26 


1.83 x 10" 26 


1.97 x 10" 26 


Main final states 


it (77.7%) 
W+W~ (8.9%) 
ZZ (6.5%) 


it (80.4%) 
W+W- (8.1%) 
ZZ (5.8%) 


W+W- (65.6%) 
Z°Z° (27.5%) 
Z°h (7.2%) 


W+W' (65.2%) 
Z°Z° (26.8%) 
Z°h (7.2%) 



Table 4: Main features of the two sets of parameters used for the comparison 



For the signal itself, we compare the results obtained with DarkSusy 5.0.4 and micrOMEGAs2.4. 
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To perform the comparison, we generate MSSM points with DarkSusy and used the result- 
ing SLHA file [39] as an input for micrOMEGAs . Each package computes the annihilation 
cross-section, branching ratios, propagation and distribution. We have chosen two sam- 
ple models in the MSSM, the input parameters at the weak scale are specified in Tab. H] 
where all dimensionful parameters are in GeV. We assume 2M\ = M 2 = M 3 /3 as would 
be obtained in a GUT model with gaugino mass universality, a common sfermion mass 
at the weak scale, mo and vanishing trilinear couplings of sleptons and of the first and 
second generations of squarks. The files used for this comparison are provided in the 
MSSM directory of micrOMEGAs(modell . slha and model2 . slha). The main features of 
the sample models are summarized in Tab.Hl The first obsevation is that the total annihi- 
lation cross section as well as the individual channel contribution can vary by up to 17%. 
This is particularly important if the main annihilation channels are into third generation 
quarks (Model 1). Indeed DarkSusy and micrOMEGAs use a different prescription for 
the running of quark Yukawa couplings. Part of the difference between the two codes is 
due to a different prescription for sin 2 9w, an On-Shell scheme in DarkSusy and MS in 
micrOMEGAs . 

To compare the 7-ray spectrum, we use an NFW profile and compute the signal 
from the Galactic center in a cone of 1° opening angle (corresponding to a solid angle of 
AO = 10~ 3 sr). Fig. [3] displays the results for both models. For this we have switched on 
the gauge boson polarisation and included the contribution from 3-body final states with 
a photon. The small difference between the two spectra is explained by the difference in 
ov. Notice that on these plots, the flux is i? 2 -corrected. 

For the one-loop processes leading to a monochromatic gamma-ray line, a detailed 
comparison between DarkSusy and SloopS was performed in [21] . The two codes agreed 
very well for the annihilation of a neutralino pair into two photons. We also find good 
agreement (< 5%) between micrOMEGAs and DarkSusy for our two test models. On 
the other hand in [2T] it was pointed out that large differences between DarkSusy and 
SloopS can occur for the 7Z one-loop process (> 30%). This is mainly because some 
diagrams were neglected in DarkSusy , these give a non negligible contribution for the 
case of a higgsino LSP. 




1 10 10 2 1 10 10 2 

E (GeV) E (GeV) 



Figure 3: 7-ray signal for Model l(left) and Model 2(right) with micrOMEGAs and 
DarkSusy 

In the case of positrons, the comparison is shown on Fig. |H The propagation pa- 
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rameters are the default parameters of DarkSusy , R = 30kpc, h = O.lkpc, L = 4kpc, 
K = 0.0826kpc 2 /Myr, 5 = 0.6, V c = lOkm/s. Note that the space diffusion coefficient K 
is defined at a reference energy of 4GeV in DarkSusy rather than lGeV in micrOMEGAs . 
To reproduce similar propagation parameters one must therefore rescale the coefficient 
K Q in micrOMEGAs. Furthermore for positrons DarkSusy assumes that 5 = below 
p = 4GeV. The difference between the two codes is small (few percent) once taking into 
account the correction due to the initial av. Larger differences are found at low energies 
and could be due to the different treatment of the positron propagation in that energy 
range. However, as shown below in section 5.3, the observed differences are much be- 
low the theoretical uncertainty induced by the propagation parameters. Good agreement 
between the two codes is also recovered for the antiproton signal, see Fig. 




1 10 10 2 1 10 10 2 



E (GeV) E (GeV) 

Figure 4: Positron signal for Model 1 (left) and Model 2 (right) with micrOMEGAs and 
DarkSusy. Here we have set 5 = 0.6, Kq = 0.03607 kpc /Myr, L = 4 kpc and Vc = 
10 km/s. 




1 10 10 2 1 10 10 2 

E (GeV) E (GeV) 



Figure 5: Antiproton signal for Model l(left) and Model 2 (right) with micrOMEGAs and 
DarkSusy. Same propagation parameters as above. Solar modulation is included with 
(j) F = 320 MV. 
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5.3 Theoretical uncertainties related to propagation 

The way the propagation of charged particles is handled in the code allows to very 
quickly estimate the propagation related uncertainty. For this we compare antiprotons 
and positrons signals using extreme values for the propagation parameters that are still 
allowed by all cosmic ray measurements. The curves labelled "Min", "Med", "Max" in 
Fig. |6] correspond to the parameters specified in Tab. [2j Recall that these parameters 
were selected to reproduce the data on the B/C ratio [301 EI]- For antiprotons the un- 
certainty exceeds one order of magnitude over the full range of energy while for positrons 
the uncertainty is large mainly at low energies. Actually high energy positrons do not 
suffer from propagation uncertainties since they are produced locally. 




Figure 6: Propagation related uncertainty for antiprotons in Model 1 (left) and for 
positrons in Model 2 (right). The propagation parameters are listed in Table 2. 



6 Conclusions 

With the new module for indirect detection presented here, micrOMEGAs is a comprehen- 
sive code for the study of the properties of weakly interacting cold dark matter candidates 
in various extensions of the standard model. The main features of this new version include 
the propagation of charged cosmic rays and some higher order processes in DM annihi- 
lation. In particular the code contains a full computation of photon radiation. These 
new features are available for any model with a weakly interacting particle that is already 
implemented in micrOMEGAs. In addition a specific code for the computation of the loop- 
induced gamma-ray line in the MSSM is included. The modular and flexible structure of 
the code makes it simple for the user to improve some specific part of the code, for exam- 
ple adding new functions that take into account the presence of dark matter clumps. The 
main observable that is not yet included in the code is the neutrino signal associated with 
dark matter capture in the Sun or in the Earth. This feature will be available in the next 



upgrade of micrOMEGAs. The code is available at http://lapth.in2p3.fr/micromegas/ 
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